/* *************************************** */
/* * Smooth - CUDA Edition               * */
/* * Authors:                            * */
/* *  André Bannwart Perina      7152757 * */
/* *  Rodrigo de Freitas Pereira 7573472 * */
/* *  Thaís Emanuele dos Santos  6453087 * */
/* *************************************** */

/* NOTE: Blocks should be at least 2x2 threads wide! */

extern "C" { 
	#include "../include/ppm_funcs.h"
}

#include <stdio.h>
#include <sys/time.h>

/* Block Dimension */
#define BLOCK_DIM 8

/* Smooth Kernel */
__global__ void cuda_smooth(byte *in, byte *out, unsigned int *m, unsigned int *n) {
	/* Block Shared Memory */
	__shared__ byte temp[BLOCK_DIM + 4][BLOCK_DIM + 4];
	/* Global Index: Thread position in global memory */
	int gIdx = ((gridDim.x * blockDim.x) * ((blockDim.y * blockIdx.y) + threadIdx.y)) + ((blockDim.x * blockIdx.x) + threadIdx.x);
	/* Thread position considering borders */
	int tX = threadIdx.x + 2, tY = threadIdx.y + 2, tDX = blockDim.x + 2, tDY = blockDim.y + 2;
	unsigned int res;

	/* Set shared memory with thread's position element */
	temp[tY][tX] = in[gIdx];

	/* Thread is at left edge, get elements from left border in the same line */
	if(0 == threadIdx.x) {
		temp[tY][0] = (0 == blockIdx.x)? 0 : in[gIdx - 2];
		temp[tY][1] = (0 == blockIdx.x)? 0 : in[gIdx - 1];
	}
	/* Thread is at right edge, get elements from right border in the same line */
	if(blockDim.x - 1 == threadIdx.x) {
		temp[tY][tDX] = (gridDim.x - 1 == blockIdx.x)? 0 : in[gIdx + 1];
		temp[tY][tDX + 1] = (gridDim.x - 1 == blockIdx.x)? 0 : in[gIdx + 2];
	}
	/* Thread is at top edge, get elements from top border in the same column */
	if(0 == threadIdx.y) {
		temp[0][tX] = (0 == blockIdx.y)? 0 : in[gIdx - (2 * *n)];
		temp[1][tX] = (0 == blockIdx.y)? 0 : in[gIdx - *n];
	}
	/* Thread is at bottom edge, get elements from bottom border in the same column */
	if(blockDim.y - 1 == threadIdx.y) {
		temp[tDY][tX] = (gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + *n];
		temp[tDY + 1][tX] = (gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + (2 * *n)];
	}
	/* Upper-left corner */
	if(0 == threadIdx.x && 0 == threadIdx.y) {
		temp[0][0] = (0 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - (2 * *n) - 2];
		temp[0][1] = (0 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - (2 * *n) - 1];
		temp[1][0] = (0 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - *n - 2];
		temp[1][1] = (0 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - *n - 1];
	}
	/* Upper-right corner */
	if(blockDim.x - 1 == threadIdx.x && 0 == threadIdx.y) {
		temp[0][tDX] = (gridDim.x - 1 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - (2 * *n) + 1];
		temp[0][tDX + 1] = (gridDim.x - 1 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - (2 * *n) + 2];
		temp[1][tDX] = (gridDim.x - 1 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - *n + 1];
		temp[1][tDX + 1] = (gridDim.x - 1 == blockIdx.x || 0 == blockIdx.y)? 0 : in[gIdx - *n + 2];
	}
	/* Lower-left corner */
	if(0 == threadIdx.x && blockDim.y - 1 == threadIdx.y) {
		temp[tDY][0] = (0 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + *n - 2];
		temp[tDY][1] = (0 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + *n - 1];
		temp[tDY + 1][0] = (0 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + (2 * *n) - 2];
		temp[tDY + 1][1] = (0 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + (2 * *n) - 1];
	}
	/* Lower-right corner */
	if(blockDim.x - 1 == threadIdx.x && blockDim.y - 1 == threadIdx.y) {
		temp[tDY][tDX] = (gridDim.x - 1 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + *n + 1];
		temp[tDY][tDX + 1] = (gridDim.x - 1 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + *n + 2];
		temp[tDY + 1][tDX] = (gridDim.x - 1 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + (2 * *n) + 1];
		temp[tDY + 1][tDX + 1] = (gridDim.x - 1 == blockIdx.x || gridDim.y - 1 == blockIdx.y)? 0 : in[gIdx + (2 * *n) + 2];
	}

	/* Wait for all threads to synchronise */
	__syncthreads();

	/* Calculate smooth */
	res =  ((temp[tY-2][tX-2] & 0xff) +
				(temp[tY-2][tX-1] & 0xff) +
				(temp[tY-2][tX] & 0xff) +
				(temp[tY-2][tX+1] & 0xff) +
				(temp[tY-2][tX+2] & 0xff) +
				(temp[tY-1][tX-2] & 0xff) +
				(temp[tY-1][tX-1] & 0xff) +
				(temp[tY-1][tX] & 0xff) +
				(temp[tY-1][tX+1] & 0xff) +
				(temp[tY-1][tX+2] & 0xff) +
				(temp[tY][tX-2] & 0xff) +
				(temp[tY][tX-1] & 0xff) +
				(temp[tY][tX] & 0xff) +
				(temp[tY][tX+1] & 0xff) +
				(temp[tY][tX+2] & 0xff) +
				(temp[tY+1][tX-2] & 0xff) +
				(temp[tY+1][tX-1] & 0xff) +
				(temp[tY+1][tX] & 0xff) +
				(temp[tY+1][tX+1] & 0xff) +
				(temp[tY+1][tX+2] & 0xff) +
				(temp[tY+2][tX-2] & 0xff) +
				(temp[tY+2][tX-1] & 0xff) +
				(temp[tY+2][tX] & 0xff) +
				(temp[tY+2][tX+1] & 0xff) +
				(temp[tY+2][tX+2] & 0xff)) / 25;

	/* Save result */
	out[gIdx] = res & 0xff;
}

int main(int argc, char *argv[]) {
	unsigned int length, width, *devLength, *devWidth;
	ppm_seg_t *r, *g, *b;
	byte *rDevChunk, *gDevChunk, *bDevChunk, *rDevOutChunk, *gDevOutChunk, *bDevOutChunk;
	struct timeval tEnd, tInit, tPEnd, tPInit;
	cudaStream_t rStream, gStream, bStream;

	gettimeofday(&tInit, NULL);

	/* Get image size */
	ppm_get_size("in.ppm", PPM_FORMAT_P3, &length, &width);

	/* Image shoud comply with block dimension */
	if(width % BLOCK_DIM || length % BLOCK_DIM) {
		fprintf(stderr, "%s: width and height should be a multiple of %d.\n", argv[0], BLOCK_DIM);

		return EXIT_FAILURE;
	}

	/* Set number of blocks and number of threads */
	dim3 nBlocks(width / BLOCK_DIM, length / BLOCK_DIM);
	dim3 nThreads(BLOCK_DIM, BLOCK_DIM);

	/* Create constant segments */
	ppm_create_const(&r);
	ppm_create_const(&g);
	ppm_create_const(&b);

	/* Get channels from file */
	ppm_get_channel_const("in.ppm", PPM_CHANNEL_R, PPM_FORMAT_P3, r);
	ppm_get_channel_const("in.ppm", PPM_CHANNEL_G, PPM_FORMAT_P3, g);
	ppm_get_channel_const("in.ppm", PPM_CHANNEL_B, PPM_FORMAT_P3, b);

	/* Allocate GPU stuff */
	cudaMalloc(&devLength, sizeof(unsigned int));
	cudaMalloc(&devWidth, sizeof(unsigned int));
	cudaMalloc(&rDevChunk, length * width);
	cudaMalloc(&rDevOutChunk, length * width);
	cudaMalloc(&gDevChunk, length * width);
	cudaMalloc(&gDevOutChunk, length * width);
	cudaMalloc(&bDevChunk, length * width);
	cudaMalloc(&bDevOutChunk, length * width);

	/* Send sizes to GPU */
	cudaMemcpy(devLength, &length, sizeof(unsigned int), cudaMemcpyHostToDevice);
	cudaMemcpy(devWidth, &width, sizeof(unsigned int), cudaMemcpyHostToDevice);

	/* Create streams */
	cudaStreamCreate(&rStream);
	cudaStreamCreate(&gStream);
	cudaStreamCreate(&bStream);

	gettimeofday(&tPInit, NULL);

	/* Copy channels to GPU */
	cudaMemcpyAsync(rDevChunk, r->d[0], length * width, cudaMemcpyHostToDevice, rStream);
	cudaMemcpyAsync(gDevChunk, g->d[0], length * width, cudaMemcpyHostToDevice, gStream);
	cudaMemcpyAsync(bDevChunk, b->d[0], length * width, cudaMemcpyHostToDevice, bStream);

	/* Apply smooth concurrently */
	cuda_smooth<<<nBlocks,nThreads,0,rStream>>>(rDevChunk, rDevOutChunk, devLength, devWidth);
	cuda_smooth<<<nBlocks,nThreads,0,gStream>>>(gDevChunk, gDevOutChunk, devLength, devWidth);
	cuda_smooth<<<nBlocks,nThreads,0,bStream>>>(bDevChunk, bDevOutChunk, devLength, devWidth);

	/* Get channels from GPU */
	cudaMemcpyAsync(r->d[0], rDevOutChunk, length * width, cudaMemcpyDeviceToHost, rStream);
	cudaMemcpyAsync(g->d[0], gDevOutChunk, length * width, cudaMemcpyDeviceToHost, gStream);
	cudaMemcpyAsync(b->d[0], bDevOutChunk, length * width, cudaMemcpyDeviceToHost, bStream);

	/* Wait for all streams to finish */
	cudaStreamSynchronize(rStream);
	cudaStreamSynchronize(gStream);
	cudaStreamSynchronize(bStream);

	gettimeofday(&tPEnd, NULL);

	printf("Elapsed time (only processing): %lu us\n", (((tPEnd.tv_sec * 1000000) + tPEnd.tv_usec) - ((tPInit.tv_sec * 1000000) + tPInit.tv_usec)));

	/* Save file */
	ppm_set_file("out.ppm", PPM_FORMAT_P3, r, g, b);

	/* Free stuff */
	cudaFree(rDevChunk);
	cudaFree(rDevOutChunk);
	cudaFree(gDevChunk);
	cudaFree(gDevOutChunk);
	cudaFree(bDevChunk);
	cudaFree(bDevOutChunk);
	cudaFree(devLength);
	cudaFree(devWidth);
	ppm_destroy_const(&r);
	ppm_destroy_const(&g);
	ppm_destroy_const(&b);
	cudaStreamDestroy(rStream);
	cudaStreamDestroy(gStream);
	cudaStreamDestroy(bStream);

	gettimeofday(&tEnd, NULL);

	printf("Total Elapsed time: %lu us\n", (((tEnd.tv_sec * 1000000) + tEnd.tv_usec) - ((tInit.tv_sec * 1000000) + tInit.tv_usec)));

	return EXIT_SUCCESS;
}